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ABSTRACT 

Analysis of the transit light curve deformed by the stellar gravity darkening allows us to photomet¬ 
rically measure both components of the spin-orbit angle ip, its sky projection A and inclination of the 
stellar spin axis z*. In this paper, we apply the method to two transiting hot Jupiter systems moni¬ 
tored with the Kepler spacecraft, Kepler-13A and HAT-P-7. For Kepler-13A, we find i k = 81° ±5° and 
ip = 60° ± 2° adopting the spectroscopic constraint A = 58?6 ± 2?0 by Johnson et al. (2014). In our 
solution, the discrepancy between the above A and that previously reported by Barnes et al. (2011) is 
solved by fitting both of the two parameters in the quadratic limb-darkening law. We also report the 
temporal variation in the orbital inclination of Kepler-13Ab, d| cosz or b|/di = (—7.0±0.4) x 10 _6 day^ 1 , 
providing further evidence for the spin-orbit precession in this system. By fitting the precession model 
to the time series of i a r b, A, and z* obtained with the gravity-darkened model, we constrain the stellar 
quadrupole moment J2 = (6.1 ± 0.3) x 10~ 5 for our new solution, which is several times smaller than 
J 2 = (1.66±0.08) x 1CU 4 obtained for the previous one. We show that the difference can be observable 
in the future evolution of A, thus providing a possibility to test our solution with follow-up obser¬ 
vations. The second target, HAT-P-7, is the first F-dwarf star analyzed with the gravity-darkening 
method. Our analysis points to a nearly pole-on configuration with ip = 101° ± 2° or 87° ± 2° and 
the gravity-darkening exponent /3 consistent with 0.25. Such an observational constraint on /3 can be 
useful for testing the theory of gravity darkening. 

Subject headings: planets and satellites: individual (Kepler-13, KOI-13, KIC 9941662) - planets and 
satellites: individual (HAT-P-7, KOI-2, KIC 10666592) - stars: rotation - tech¬ 
niques: photometric 


1. INTRODUCTION 

Spin-orbit angle or the stellar obliquity, ip , the angle 
between the stellar spin axis and the orbital axis of its 
planet, serves as a unique probe of the dynamical history 
of planetary systems. Especially, its connection with the 
hot-Jupiter migration has been extensively studied (e.g., 
Fabrycky & Winn 2009), but the relationship between 
the observed samples and the migration process is not 
straightforward for various reasons. First of all, the ini¬ 
tial distribution of the spin-orbit angles is not known. 
Some studies do suggest that the protoplanetary disk 
may have already been misaligned with the stellar equa¬ 
tor due to the chaotic gas accretion (e.g., Bate et al. 
2010; Fielding et al. 2014) or the magnetic star-planet 
interaction (e.g., Lai et al. 2011). In these cases, the 
spin-orbit misalignment is primordial, rather than due 
to the migration. Even after the disk dissipation or the 
completion of migration, spin-orbit angle can evolve due 
to the gravitational perturbation from the companion 
(e.g., Storch et al. 2014; Li et al. 2014). As suggested 
by the observed correlation between the spin-orbit mis¬ 
alignment and stellar effective temperature (Winn et al. 
2010; Albrecht et al. 2012), spin-orbit angle may also be 
affected by the tidal star-planet interaction (e.g., Xue 
et al. 2014), whose mechanism is not well understood. 
To partially resolve these issues, it is beneficial to mea¬ 
sure spin-orbit angles for systems with various host-star 
and orbital properties. For instance, planets on distant 
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orbits or around hot/young stars are valuable targets be¬ 
cause we expect that tides have not significantly affected 
the primordial spin-orbit configuration. 

This paper focuses on a relatively new method for 
the spin-orbit angle determination in transiting systems, 
which utilizes the gravity darkening of the host star ow¬ 
ing to its rapid rotation (Barnes 2009). Stellar rotation 
makes the effective surface gravity at the stellar equa¬ 
tor smaller than that at the pole by a fractional or¬ 
der of 7 = fllR3/2GM k ~ (P br /P ro t) 2 , where H*, P*, 
M*., .Pbr, and P ro t are angular rotation frequency, radius, 
mass, break-up rotation period, and rotation period of 
the star, respectively. According to von Zeipel’s theorem 
(von Zeipel 1924), this results in the inhomogeneity of the 

stellar surface brightness through the relation T e g oc g^ s . 
Here, P e ff and g e s are the effective temperature and sur¬ 
face gravity at each point on the stellar surface, and 
gravity-darkening exponent /3 characterizes the strength 
of the gravity darkening, which is theoretically 0.25 for a 
barotropic star with a radiative envelope. When a planet 
transits a star with such an inhomogeneous and generally 
non-axisymmetric brightness distribution, an anomaly of 
0(76) appears in the light curve, where 8 is the transit 
depth. Since the shape of the anomaly depends on the 
position of the stellar pole relative to the planetary or¬ 
bit, the obliquity of the stellar spin ip can be estimated 
with the light-curve model taking into account the effect 
of gravity darkening. 

Indeed, this “gravity-darkening method” has many 
unique aspects. So far, it is the only known method 
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that simultaneously constrains both components of ip, the 
sky-projected spin-orbit angle A and stellar inclination 
i* (c.f., Equation 2 and Figure 1). Moreover, obliquity 
analysis is possible essentially with the photometric data 
alone, and its application is not necessarily limited to 
short-period planets, as far as the transit is observed with 
sufficient signal-to-noise ratio (Zhou & Huang 2013). It 
is also interesting to note that the method is (only) ap¬ 
plicable to fast-rotating (i.e., young or hot) stars, for 
which anomalies of larger amplitudes result. Since rapid 
rotators are not suitable for the precise spectroscopic 
velocimetry because of their broad spectral lines, this 
method is complementary to the conventional spin-orbit 
angle measurement using the Rossiter-McLaughlin (RM) 
effect. All these properties make the method suitable for 
sampling stars for which tidal effect is not so significant 
that the primordial information is expected to be well 
preserved in the current spin-orbit configuration. 

Although the gravity-darkening method is valuable in 
many aspects, the procedure for obtaining ip may not be 
fully established. In a representative example of its ap¬ 
plication, Kepler-13A, the constraint from the gravity- 
darkening method (Barnes et al. 2011, hereafter Bll) 
is known to be in disagreement with the later spectro¬ 
scopic measurement of A with the Doppler tomography 
(Johnson et al. 2014). In addition, inconsistent results 
arise even within the gravity-darkening analyses, depend¬ 
ing on the choice of the limb-darkening coefficients or /3 
(Zhou & Huang 2013; Ahlers et al. 2014). For these rea¬ 
sons, it is worth revisiting the reliability and limitation 
of this method more carefully, in order for this unique 
method to be applied to more systems in future and pro¬ 
vide credible results. 

In this paper, we reanalyze a well-known example of 
the gravity-darkened transit of Kepler-13Ab, with more 
data than used in the previous analysis by Bll. We 
investigate the systematic effects in the spin-orbit an¬ 
gle determination, and propose a joint solution that may 
solve the discrepancy with the Doppler tomography mea¬ 
surement (Section 3). We will also see that the spin-orbit 
precession in this system can be used to test the valid¬ 
ity of our solution, as well as to determine the stellar 
quadrupole moment </ 2 (Section 4). 

In addition, we apply the gravity-darkening method 
for the first time to an F-type dwarf star, HAT-P-7, 
where the anomaly in the transit light curve has been 
reported in several studies (e.g., Esteves et al. 2013; Van 
Eylen et al. 2013; Esteves et al. 2014; Benomar et al. 
2014). While the RM measurements (Winn et al. 2009a; 
Narita et al. 2009; Albrecht et al. 2012 ) have established 
that A > 90°, suggesting a retrograde orbit, the follow¬ 
ing asteroseismic inferences (Benomar et al. 2014; Lund 
et al. 2014) have revealed that a pole-on orbit is actually 
favored. In Section 5, we show that a similar conclu¬ 
sion is also obtained from the gravity-darkening method 
and discuss the consistency of our result with other con¬ 
straints on the host-star properties. 

2 . METHOD 
2.1. Model 

We basically follow Barnes (2009) in modeling the 
gravity-darkened transit light curve. The model includes 
the following 14 parameters, which are listed as “fitting 


parameters” in Tables 1 and 3: 

1 . mean stellar density, p* = 3 -M*/ 47 ri?*, which corre¬ 
sponds to the semi-major axis scaled by the stellar 
equatorial radius, a/i ?* 1 

2 . limb-darkening coefficient for the quadratic law, 
Cl = Ml + m 2 , 

3. limb-darkening coefficient for the quadratic law, 
c 2 = Mi - m 2 , 

4. time of the inferior conjunction, t c , 

5. orbital period, P, 

6 . cosine of orbital inclination, cosf or bi 

7. planetary radius normalized to the stellar equato¬ 
rial radius, i? p /i?* 

8 . normalization of the out-of-transit flux, Fq 

9. stellar mass, M*, 

10 . stellar rotation frequency, / ro t 

11 . stellar effective temperature at the pole, T* iPO i e 

12 . gravity-darkening exponent, /3, 

13. stellar inclination, i* 

14. sky-projected spin-orbit angle, A. 

The first eight parameters are common with the light- 
curve model without gravity darkening. We assume cir¬ 
cular orbits for the two targets because the orbital eccen¬ 
tricities are constrained to be very small, if any, from the 
occultation light curves (Shporer et al. 2014; Benomar 
et al. 2014). 

In the gravity-darkened model by Barnes (2009), the 
shape of the star is approximated by the spheroid 
with the oblateness 7 = fi 2 i?3/2GM* = 37 r/ 2 ot / 2 Gp*. 
The surface brightness at each point is modeled as 
the blackbody emission of the temperature T* = 

T 1 *,pole ( 5 eff/Pefi',poie) /3 , where g e s/g e s, pole is the effective 
surface gravity normalized by its value at the stellar pole. 
The surface gravity at point r on the stellar surface is 
calculated by g e s = —GM+r~ 2 r + 4tt 2 fp ot r±f±. Here r 
and r are the norm and unit vector of the radius vector 
r, respectively. Similarly, r± and r± are those of r±, 
the projection of r onto the stellar equatorial plane. The 
Planck function B\(T+) at each point is convolved with 
the “high-resolution” Kepler response function 2 using the 
table of the wavelength- and temperature-dependent fac¬ 
tor calculated prior to the fitting. The convolved flux is 
then multiplied by the limb-darkening function 

I{h) = 1 - w(i - h) ~ m 2 (1 - p) 2 , (1) 

with fi being the cosine of the angle between —g e ff and 
our line of sight, 3 and integrated over the visible surface 

1 In this paper, i?* denotes the equatorial radius of the star. 

2 http://keplergo.arc.nasa.gov/CalibrationResponse.shtml 

3 Although this vector —g e ff is not exactly parallel to the surface 
normal of the spheroid we assume, the difference is 0( j 2 ) and thus 
negligible. 
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Fig. 1.— Definitions of i 0 rtn *★> A, and 0 in this paper. The 
orbital inclination, i Q r b, is the angle between the planetary orbital 
axis (blue arrow) and the observer’s line of sight. In a transiting 
system, i or b is usually very close to 7r/2 and hence the orbital 
axis almost coincides with its projection onto the plane of the sky. 
Inclination of the stellar spin axis, ?*, is similarly defined as the 
angle between the stellar spin axis (red arrow) and the line of sight. 
The angle between the two axes (red and blue ones), 0, is the spin- 
orbit angle or the stellar obliquity. Its sky projection, A, denotes 
the angle between the sky projections of the same two axes. 


of the star to give the total flux. We fix T* iPO i e at the ob¬ 
served effective temperature assuming that the difference 
between T* iPO i e and the disk-integrated effective temper¬ 
ature is small. Note that the gravity-darkened transit 
light curve gives p* alone and can not constrain M* and 
R+ separately, as is the case for the transit without grav¬ 
ity darkening. 

The configuration of the planetary orbit and stellar 
spin is specified by three angles, Arbi i *, and A, which 
are defined in Figure 1 (see also figure 1 of Benomar et al. 
2014). The orbital and stellar inclinations, Arb and A> 
are measured from the line of sight and defined to be in 
the range [0,7r]. The sky-projected spin-orbit angle, A, 
is the angle between the sky-projected stellar spin and 
planetary orbital axes. It is measured from the former to 
the latter counterclockwise in the sky plane, and is in the 
range [0, 27 t] . With these definitions, the true spin-orbit 
angle, or the stellar obliquity, ip, is given by equation 1 
of Benomar et al. (2014): 

cos ip = cos A cos Arb + sin i * sin Arb cos A. (2) 


This, in principle, allows us to break the degeneracy be¬ 
tween M* and I?*, enabling the determination of the ab¬ 
solute dimension of the system. Nevertheless, the con¬ 
straint on M„ is usually weak as discussed in Bll, and 
so we fix M* at the observed value. 

2.2. Data Processing 

We detrend and normalize the transit light curves of 
each target along with the consistent determination of 
the transit times and transit parameters. We first nor¬ 
malize the light curve of each quarter using its median, 
and then iterate the following two steps until the result¬ 
ing transit times t r and transit parameters converge (typ¬ 
ically 10-20 times): 

1. Light curve around each transit (±0.2 days for 
Kepler-13A and ±0.15 days for HAT-P-7) is mod¬ 
eled as the product of a quadratic polynomial 4 
a-o ± a\(t — t c ) ± 02 (t — t c ) 2 (t: time) and the an¬ 
alytic transit light-curve model by Mandel & Agol 
(2002). We use the Levenberg-Markwardt (LM) 
method (Markwardt 2009) to fit ao, a\, a%, and 
t c iteratively removing 5a outliers, while the other 
parameters are fixed. The filtered data are then di¬ 
vided by the best-fit polynomial to give a normal¬ 
ized and detrended transit light curve. We discard 
the transits with data gaps of more than 50%. 

2. Using the set of t c obtained in the first step, we 
calculate the mean orbital period P and transit 
epoch to by linear fit and use them to phase-fold 
the normalized and detrended transits. The phase- 
folded light curve is averaged into one-minute bin 
and then fitted with the Mandel & Agol (2002) 
model using an LM algorithm. We fit ci, C 2 , p*, 
cos forb, Rp/R*, and F 0 , whose best-fit values are 
used in the step 1 of the next iteration. In this 
step, the orbital period P is fixed to be the value 
obtained from the linear fit and the central time of 
the phase-folded transit is fixed to be zero. 

In the following analysis, we use the one-minute 
binned, phase-folded light curve obtained in the second 
step of the final iteration. For each bin, the flux value 
is given by its mean and the error is estimated as the 
standard deviation within the bin divided by the square 
root of the number of data points. 


Throughout the paper, we restrict A to be in the range 
[0,7 t/ 2] making use of the intrinsic symmetry with re¬ 
spect to the sky plane. We do not lose any physi¬ 
cal information of the system with this choice because 
any of the relative star-planet configurations with A in 
[7 t/2,7t] is the same as one of those with A in [0,7r/2]. 
In other words, the configurations (U, Wb, A) and (n — 
A, 7r — Arb, —A) are equivalent. This transformation cor¬ 
responds to looking at the system from the other side of 
the plane of the sky. 

In the following, we also adopt the constraint on the 
stellar line-of-sight rotational velocity v sin A from spec¬ 
troscopy, which is related to the above model parameters 

by 

. . _ , /3 M* 

v sin = 27r/ rot -- 

\4irp* 



2.3. Fitting Procedure 

In fitting the observed light curves, the likelihood C of 
the model is computed by C oc exp(—% 2 /2), where 

- \ 2 
Pj Pmodelj \ 

S Pj ) 

In the first term, /), /model, u and cq are the observed 
value, modeled value, and error of the ith flux data. The 
second term is introduced to take into account the con¬ 
straints from other observations on some (functions) of 
the model parameters pj. In the following analysis, p is 

4 Use of the quadratic polynomial helps the better removal of 
flux variation not due to the transit, i.e., planetary light, ellipsoidal 
variation, and Doppler beaming. 
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read to be v sin t* and, in some cases, A. 5 For each pj , we 
assume a Gaussian constraint of the form pj ± Spj and 
the value obtained from the model is denoted by p m odei, j ■ 

The maximum likelihood solution is found by minimiz¬ 
ing Equation (4) with the LM method using the empfit 
package (Markwardt 2009). Since the complex depen¬ 
dence of x 2 011 i* and A is expected, we repeat the fit¬ 
ting procedure from the initial i* in [0,90°] and A in 
[—180°, 180°] at 10° intervals. Initial values of the other 
parameters are chosen close to the best-fit values ob¬ 
tained from the model without gravity darkening. We 
also try both positive and negative cosi or b as an initial 
value to search the whole domain of i 0 rb, which is now 
[0°, 180°]. 

3. TRANSIT ANALYSIS OF KEPLER-13Ab 

In this section, we report the analysis of the gravity- 
darkened transit of Kepler-13Ab. We first analyze the 
whole available short-cadence (SC) data from Q2, 3, 
and 7-17 using the same stellar parameters as in Bll 
to test the validity of our method (Section 3.1). Moti¬ 
vated by the recently reported disagreement with A from 
the Doppler tomography, we also investigate the possible 
systematics in the spin-orbit determination arising from 
the choice of stellar parameters. We show that the dis¬ 
crepancy can be absorbed by adjusting the value of C2 
and present a joint solution that is compatible with all 
of the observations made so far. 

3.1. Reproducing the Results by Bll 

In this subsection, we analyze the short-cadence (SC), 
Pre-search Data Conditioned Simple Aperture Photom¬ 
etry (PDCSAP) fluxes from Q2, 3, and 7 17. Given the 
clear transit duration variation (TDV) reported by Szabo 
et al. (2012) and Szabo et al. (2014), we separately an¬ 
alyze the transits from each quarter, rather than folding 
all the available data. Since we do not detect significant 
temporal variations in the parameters other than cos i or b 
(see Section 4), we report the mean and standard devia¬ 
tion of the best-fit values from the above 13 quarters for 
each parameter. 

First, we use the same stellar parameters as in Bll 
and obtain the results in the second column of Table 1. 
Namely, we subtract a constant value F c = 0.45 from the 
normalized flux to remove the flux contamination from 
the companion star, and impose the constraint usini* = 
GSilOkms^ 1 based on Szabo et al. (2011). We fix AT* = 
1.83M 0 and T*, po i e = 8848K from Borucki et al. (2011), 
and C2 = 0. In Figure 2, the best-fit model is overplotted 
with the data for Q2, which is to be compared with figure 
2 of Bll. 

Basically, we find a very good agreement with the re¬ 
sult by Bll using about 12 times more data. Although 
the values of cosi or b> **, and A we report here appear 
different from those in Bll, that is simply because we 
choose to be in the range [0,7r/2]. This is physically 
the same configuration as theirs and corresponds to the 
top-left situation in figure 3 of Bll. That is, A in our 
solutions with cosf or b < 0 should be read as —A in the 
conventional definition, because A is usually defined for 

5 Only in Section 4.1, p*, ci, C 2 , and / ro t are also in¬ 

cluded. 


the orbit with cosi or b > 0 (see also the discussion after 
Equation 2). 

In addition to the solution in Table 1, we also find a 
retrograde solution with A > 90° as noted in Bll. Here 
we do not discuss this solution, however, because the 
Doppler tomography observation has already excluded 
the retrograde orbit with high significance (Johnson et al. 
2014). 

3.2. Systematics due to Stellar Parameters 

Although we find consistent values of A and i* as 
obtained by Bll, those of A significantly differ from 
A = 58?6 ± 2?0, the value obtained from the Doppler to¬ 
mography (Johnson et al. 2014). Motivated by this dis¬ 
crepancy, we investigate the possible origins of system¬ 
atics in the spin-orbit angle determination with gravity 
darkening in this subsection. 

First, we examine the systematics due to the choice of 
AT*, usini*, T* jPO i e , and F c , which are the stellar prop¬ 
erties not derived from the light curve modeling. 6 We 
perform the same analysis as in Section 3.1, but adopt¬ 
ing the following parameters from the most recent photo¬ 
metric and spectroscopic study by Shporer et al. (2014, 
hereafter S14): usini* = 78 ± 15kms _1 , AL* = 1.72 ATq, 
7*,pole = 7650 K, and F c = 0.47726. The correspond¬ 
ing results are shown in the third column of Table 1. 
We find that i* and A can differ by as large as 10° due 
to the choice of the above parameters, but the differ¬ 
ence is not so large as to explain the disagreement with 
the Doppler tomography. The main difference from the 
Bll case with this new set of parameters is the differ¬ 
ent constraint on / rot sin**, which is proportional to the 
combination (p*/M*) 1 / 3 usini* (c.f., Equation 3). With 
smaller AT* and larger usin**, the stellar rotation rate 
slightly higher than the Bll case is favored. We find that 
the difference in T* jPO i e is less important compared to the 
above effect. We also find that larger F c yields larger 
R v /R„ , which makes the impact parameter or | cos i or b | 
smaller to give the same ingress/egress duration. 

Next, we allow C2 = u\ — U 2 to be free, and find that 
the resulting spin-orbit angle is very sensitive to this pa¬ 
rameter. When C2 is floated, the constraints on i* and 
A become much weaker than the C2 = 0 case, as shown 
in the fourth and fifth columns of Table 1. The strong 
dependence on C2 is illustrated in Figure 3, which shows 
that A and i* vary by several tens of degrees depending on 
C2- In fact, the result indicates that the gravity-darkened 
light curve is actually compatible with the Doppler to¬ 
mography solution if we choose C2 ~ 0.25; such a solution 
will be discussed in Section 3.3. 

3.3. Joint Solution 

In Section 3.2, we found that the gravity-darkened 
light curve is compatible with the value of A estimated 
from the Doppler tomography if C2 ~ 0.25. Thus we 
repeat the analysis treating C2 as a free parameter for 
both stellar parameters by Bll and S14, but this time 
imposing additional constraint A = 58? 6 ± 2?0 from the 
Doppler tomography. The results are summarized in 

6 We do not examine the dependence on j3 here because Bll 
have already shown that a different choice of /3 = 0.19, suggested 
by the interferometric observation of Altair (Monnier et al. 2007), 
does not change the result significantly. 
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Fig. 2.— Fitting the gravity-darkened model to the Q2 transit of Kepler-13Ab. (Middle) Black dots are the phase-folded and binned 
fluxes from Q2. The thick red line shows the best-fit gravity-darkened model, while the thin blue line is the best-fit model without gravity 
darkening. (Bottom) Black dots are the residual of the best-fit gravity-darkened model. Gray open circles are those for the joint solution, 
where C 2 is fitted with the constraint A = 58?6=t2?0 from the Doppler tomography. (Top) Black dots are the residuals of the best-fit model 
without gravity darkening. Thick red line is the difference between the best-fit model with gravity darkening and that without gravity 
darkening. Dashed red line shows the same result for the joint solution. The difference between the two gravity-darkened solutions is only 
barely visible just after the ingress and before the egress. 


the last two columns in Table 1. The resulting value of 
ii, = 81° ± 5° indicates that the star is close to equator- 
on, and i)j = 60° ± 2° is slightly larger than the previous 
estimate. In terms of x^ in , these solutions equally well 
reproduce the transit anomaly as the solutions discussed 
so far, and still they are consistent with the Doppler to¬ 
mography result. Moreover, we obtain a slightly longer 
P rot , which better agrees with P rot = 25.43 ± 0.05 h es¬ 
timated by Szabo et al. (2012) and Szabo et al. (2014) 
than the solution with the gravity darkening alone. For 
these reasons, the joint solution is most favored from the 
current observations. 

We note, however, that the likelihood for the joint solu¬ 
tion is not so high as to statistically justify the introduc¬ 
tion of the additional free parameter C 2 . Furthermore, 
the plausibility of the value of C 2 in our joint solution is 
theoretically unclear. We obtain the theoretical values of 
Ci jt h — 0.6 and C 2 ,th — 0.0 from the table of Sing (2010) 
if we adopt the effective temperature and surface gravity 
by S14. Hence the value of C 2 from our joint solution 
is discrepant from C 2 ,th) they could even have opposite 


signs depending on the stellar parameters. Nevertheless, 
it is also true that theoretical values often disagree with 
the observed ones (e.g., Southworth 2008); in fact, Ci in 
the light-curve solution with C 2 = 0 is also different from 
Cyth- Therefore, we do not consider the possible devia¬ 
tions from the theoretical values crucial, and regard it as 
an open question. 7 An alternative approach to indepen¬ 
dently assess the validity of our solution is discussed in 
the next section. 

4. SPIN-ORBIT PRECESSION IN THE 
KEPLER-13A SYSTEM 

The shape of Kepler-13Ab’s transit is known to ex¬ 
hibit a long-term variation, which is likely due to the 
spin-orbit precession induced by the quadrupole mo¬ 
ment of the rapidly rotating host star (Szabo et al. 
2012, 2014). Indeed, we find the monotonic decrease in 
| cos ?' 0 rb | from the quarter-by-quarter analysis in Section 

7 For reference, we find C 2 = 0.1-0.2 if we adopt the model 
without gravity darkening (Mandel & Agol 2002), which suggests 
that the choice of C 2 = 0 is not indispensable. 
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Fig. 3.— Constraints on (A. ?*) from the gravity-darkened transit of Kepler-13Ab for the different choices of C 2 - In this illustration, data 
from Q2 are used and stellar parameters from Bll are adopted. The solid, dashed, and dotted contours respectively show la, 2cr, and 3a 
confidence regions for (A, i*) obtained from 200000 Markov Chain Monte Carlo (MCMC) samples for three fixed values of C 2 (0, 0.12, and 
0.25). The shaded areas bounded by the vertical solid, dashed, and dotted lines respectively denote la, 2a, and 3a confidence regions for A 
obtained from the Doppler tomography (Johnson et al. 2014). The sign of A is opposite to their quoted value because we are now dealing 
with the solution with cosz or b < 0 (i.e., it/2 < i or b < 7r); see also the discussion in the third paragraph of Section 3.1. 


3; the const ant-value model is rejected at the p -value of 
0.5% for this parameter using a simple x 2 test. On the 
other hand, the other model parameters are found to be 
consistent with the constant value using the same crite¬ 
rion. Therefore, our analysis confirms that the observed 
TDVs are actually due to the variation in cos/ or b, 8 fur¬ 
ther supporting the precession scenario with the more 
realistic model of the asymmetric transit light curve. 

In this section, we further examine this scenario with 
the gravity-darkened transit model. Unlike the above 
previous studies (Szabo et al. 2012, 2014) that focused 
on i or b, the gravity-darkened model allows us to addition¬ 
ally study the (non-)variations in the other two angles, A 
and «*, which should also be induced if the system is pro¬ 
cessing. 9 By fitting the analytic precession model to the 
time series of cosi 0 rbj A, and obtained from the light 
curves, we constrain the stellar quadrupole moment J 2 
and its moment of inertia coefficient C. On the basis of 
these constraints, we predict the future evolution of the 
system configuration and argue that the follow-up obser¬ 
vations of such a long-term modulation can distinguish 

8 Note that, in Szabo et al. (2012), the degeneracy between a/it* 
(or p*) and cos z or |, was not solved. 

9 If either of the angular momenta of the stellar spin or the 
orbital motion dominates, i or b or i* is almost constant. In the 
Kepler-13A system, the two angular momenta have comparable 
magnitudes and so all three angles modulate due to the preces¬ 
sion. A similar case, the PTFO 8-8695 system, has been studied 
by Barnes et al. (2013) and S. Kamiaka et al. (2015, in prepara¬ 
tion). 


the light-curve and joint solutions discussed in Section 
3. In the following, we mainly discuss the results ob¬ 
tained with the stellar parameters from S14, though the 
conclusions remain the same for the Bll parameters. 

4.1. Model parameters from each transit 

To examine the temporal variations in cos* 0 rb) **, and 
A, we fit individual transit light curves, rather than the 
phase-folded ones, for these parameters. We use the same 
two models (“light-curve solution” with C 2 = 0 and “joint 
solution” with C 2 fitted) as discussed in Section 3. In or¬ 
der not to underestimate the errors in the three angles, 
we fit all the other model parameters, p*, Ci, C 2 (for 
the joint model), t c , R p /R +, f ro t, and Fq as well, which 
should not vary temporally in our model. Using the best 
values in Table 1, we impose the constraints on these pa¬ 
rameters except for t c and To , through the second term of 
Equation (4). In fitting much noisier individual transits, 
this prescription assures that the parameters converge to 
the values consistent with those from the phase-folded 
light curves, while preserving their differences from tran¬ 
sit to transit. We also discard the transits for which 
fit does not converge due to the data gaps and/or flare¬ 
like brightening features sometimes found in the light 
curves. The resulting sequences of the transit parame¬ 
ters are plotted in Figure 4. 

As mentioned above, we again find the clear linear 
trend in cosi or b from individual transits. We fit the lin¬ 
ear model to the time series of cos/ or b using a Markov 
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TABLE 1 

Results for the transit of Kepler- 13Ab 



light-curve solution (02 = 0) 

light-curve solution (c 2 fitted) 

joint solution (c 2 fitted) 

Ref. for F c , asini*, M*, T* !po i e 

Bll 

S14 

Bll 

S14 

Bll 

S14 

(Assumed Flux Contamination) 
F c 

0.45 

0.47726 

0.45 

0.47726 

0.45 

0.47726 

( Constraints) 
t; sin (kms -1 ) 

65 ± 10 

78 ± 15 

65 ± 10 

78 ± 15 

65 ± 10 

78 ± 15 

A (deg) 





58.6 ±2.0* 

58.6 ±2.0* 

{Fitting Parameters) 

M* (M 0 ) 

1.83 (fixed) 

1.72 (fixed) 

1.83 (fixed) 

1.72 (fixed) 

1.83 (fixed) 

1.72 (fixed) 

T*, po le (K) 

8848 (fixed) 

7650 (fixed) 

8848 (fixed) 

7650 (fixed) 

8848 (fixed) 

7650 (fixed) 

p* (gem 3 ) 

0.533 ± 0.005 

0.550 ± 0.005 

0.530 ± 0.005 

0.547 ± 0.006 

0.525 ± 0.005 

0.538 ± 0.006 

Cl 

0.496 ± 0.008 

0.493 ± 0.008 

0.50 ±0.04 

0.51 ±0.02 

0.523 ± 0.005 

0.528 ±0.006 

C2 

0 (fixed) 

0 (fixed) 

0.02 ±0.28 

0.12 ±0.13 

0.20 ±0.02 

0.26 ±0.04 

t c (10 5 day)** 

—3 ± 1 

—3 ± 1 

-8 ±7 

-10 ± 7 

-9 ±4 

-11 ± 5 

P (day) 






-0.066 ± 0.004 

-0.057 ± 0.004 

-0.066 ±0.003 

-0.055 ± 0.004 

-0.064 ± 0.004 

-0.054 ± 0.004 

COS i Q rb 

Rp/R * 

0.0845 ± 0.0002 

0.0864 ± 0.0002 

0.0845 ± 0.0003 

0.0865 ± 0.0002 

0.0846 ± 0.0002 

0.0864 ± 0.0003 

F 0 







/rot (/iHz) 

12.9 ±0.4 

14.5 ± 0.6 

12.8 ±4.2 

12.9 ± 1.8 

10.2 ± 0.6 

11.6 ± 1.0 

i * (deg) 

47 ±3 

56 ±3 

60 ±20 

71 ± 16 

73 ±5 

81 ±5 

A (deg) 

-20.3 ± 1.3 

-13.9 ±1.3 

-33 ± 13 

-30 ± 12 

-58.4 ±2.0*** 

-58.5 ±2.0*** 

P 

0.25 (fixed) 

0.25 (fixed) 

0.25 (fixed) 

0.25 (fixed) 

0.25 (fixed) 

0.25 (fixed) 

( Derived Parameters) 

Prot (hr) 

21.5 ±0.7 

19.1 ± 0.8 

23 ±5 

22 ±3 

27 ±2 

24 ±2 

ip (deg) 

50 ±3 

40 ±3 

52 ±9 

42 ±6 

61 ± 2*** 

60 ±2*** 

impact parameter 

0.29 ±0.02 

0.26 ±0.02 

0.29 ±0.01 

0.25 ±0.02 

0.28 ± 0.02 

0.24 ± 0.02 

stellar oblateness 

0.022 ±0.001 

0.027 ±0.002 

0.02 ±0.02 

0.022 ± 0.006 

0.014 ±0.002 

0.018 ±0.003 

Xmin/ dof 

250/241 

249/241 

247/240 

245/240 

248/241 

246/241 


Note. — The quoted best-fit values and uncertainties are averages and standard deviations of the best-fit values obtained from 
13 quarters analyzed here. The value of Xmin is also the average of the minimum y 2 among quarters. 


* This should be read as — 58/6 ±2/0 for the solution with cos i or b < 0 discussed in this table. 

** Measured from the transit epoch to(BJD — 2454833) = 120.566 ±0.001 obtained with the transit model without gravity darkening. 

*** For A in the joint solutions, we quote the uncertainty in the constraint from the Doppler tomography. This is because the value of 
A is completely determined by this constraint and its standard deviation (several 0.1°) is not a good measure of the actual uncertainty. 
Accordingly, the quoted uncertainty in ip is also increased by taking a quadratic sum of its standard deviation and the additional 
scatter coming from the uncertainty of 2° in A. 


Chain Monte Carlo (MCMC) algorithm and obtain the 
rates of change in the upper part of Table 2. Here we only 
report the slopes for absolute values of cos* or b because 
its actual sign depends on the sign of cost*, which can 
never be determined with the current observations (we 
arbitrarily choose cos i* > 0 in this paper, as discussed af¬ 
ter Equation 2). Comparing the light-curve solution and 
joint solution, we find that the rate of | cost or b| change is 
insensitive to A or C 2 because | cos * or b | is mainly deter¬ 
mined from the transit duration. With a/i?* calculated 
from p*, our value for d| cosf or b|/dt is found to be con¬ 
sistent with db/dt = (—4.4± 1.2) x 10~ 5 day -1 by Szabo 
et al. (2012), but our constraint is several times better. 

Figure 4 also shows the abrupt systematic changes in 
i? p /i?*. These changes occur exactly in phase with the 
border of different quarters indicated with different colors 
(black and gray). For this reason, they are unlikely to be 
of physical origin, but are probably due to the seasonal 
transit depth variations similar to those reported by Van 
Eylen et al. (2013) for HAT-P-7. In addition, some of the 
parameters (most notably p* and / ro t) show the long¬ 
term modulation of the period ~ 400 days. Origins of 
these systematics are beyond the scope of this paper, 
and they are just treated as the additional scatter in the 
data. 


4.2. Fit to the observed angles and future prediction 

Among the observed time series of transit parameters 
in Figure 4, those of cos* 0 rb> A, and i k are fitted us¬ 
ing an MCMC algorithm to observationally constrain J 2 
and C. We utilize the same analytic precession model 
as in Barnes et al. (2013), which constitutes an analytic 
solution of the secular equations of motion derived by 
Boue & Laskar (2009). In this model, the orbital and 
spin angular momenta precess around the total angular 
momentum at the same angular rate given by 


D = D p 



+ sin 2 ip, 


( 5 ) 


where fi p is the precession rate of the orbital angular 
momentum around the stellar spin, and explicitly given 

by 

^ = 4/^-^) 2 ( 6 ) 

with J 2 being the stellar quadrupole moment. In the 
Kepler-13A system, the spin angular momentum, S , is 
comparable to the orbital one, L, owing to the small 
semi-major axis and rapid stellar rotation. As a conse- 
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BJD-2454833 (day) BJD-2454833 (day) 


Fig. 4.— Best-fit model parameters from each transit. The left panels are the results for the light-curve solution with C 2 = 0, while the 
right ones are for the joint solution. Errors are from the outputs of the cmpf it package. Parameters from even quarters (2, 8, 10, 12, 14, 
and 16) are shown in black, while those from odd quarters (3, 7, 9, 11, 13, 15, and 17) are in gray. For the times of inferior conjunctions, 
t c , the residuals of the linear fit (i.e., TTVs) are plotted for clarity. Solid lines in cosi or b panels are the best-fit linear models. 


quence, Cl also depends on the ratio of the two, 
L 1M P 1 /o\ 2 

S ~ \R~J ’ 


(7) 


where C is the moment of inertia coefficient of the host 
star. Thus, the independent model parameters are p*, 
/rot) ^ 2 , C, P, Mp/M*, and three angles cos* or b, A, i* at 
some epoch (here taken to be BJD = 2454833 + 800). We 
do not relate Ji to the other parameters like the stellar 
oblateness as done in Barnes et al. (2013). 

To realistically evaluate the credible intervals of J 2 


and C by marginalization, uncertainties in p*, / ro t, P, 
and Afp/Af* should also be taken into account. How¬ 
ever, these parameters are not well determined from the 
data of cosiorb, A, and i*. Thus, they are floated with 
the following Gaussian priors. The first three are as¬ 
signed the same central values and widths as in Table 
1. For the mass ratio, we take the mean and standard 
deviation of the results reported by S14, Esteves et al. 
(2014), and Faigler & Mazeh (2014), which come from 
the amplitudes of the ellipsoidal variation and Doppler 
beaming. We also impose the Gaussian prior on C cen- 
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TABLE 2 

Results of the precession model fit to cosi or b, i*, and A from each transit 



light-curve 

solution (02 = 0 ) 

joint solution (c 2 fitted) 

Ref. for F c , asm**, M+, T* !po i e 

Bll 

S14 

Bll 

S14 

(Linear fit to cosi or b) 

1 COSlorbf 

0.0668 ± 0.0001 

0.0581 ± 0.0001 

0.0658 ± 0.0001 

0.0560 ± 0.0002 

d ' c°d* 0rb ! (day -1 ) 

(-5.9 ±0.3) x 10 - 

6 (— 6.8 ± 0.3) x 10 -6 

(-6.0 ±0.3) x 10 -6 

(-7.0 ±0.4) x 10 -6 

(Precession model fit to cosi or b, 

i*, and A) 




P*, /rot, P 


Same as Table 1 (priors = posteriors) 


COS i 0 rb 

-0.0668 ± 0.0001 

-0.0581 ± 0.0001 

-0.0658 ± 0.0001 

-0.0560 ± 0.0002 

i* (deg)* 

44.7 ±0.3 

54.2 ± 0.3 

72.8 ±0.3 

81.8 ± 0.2 

A (deg)* 

- 20.1 ± 0.2 

-13.9 ±0.1 

-58.65 ± 0.09 

-58.62 ±0.09 

M p /M*** 

(3.4 ±0.8) x 10 -a 

1 ( 2.8 ± 0 . 8 ) x 10 -3 

(4.1 ±0.8) x 10 -3 

(4.0 ±0.8) x 10 -3 

c*** 

0.09 ± 0.02 

0.10 ± 0.02 

0.08 ± 0.02 

0.08 ± 0.02 

J 2 

(1.44 ±0.07) x IQ - 

4 ( 1.66 ± 0.08) x 10 -4 

(5.6 ±0.3) x 10 -5 

(6.1 ± 0.3) x 10 -5 

(Derived from the precession model) 




Precession period (yr) 

(5.7 ±0.4) x 10 2 

(4.3 ±0.3) x 10 2 

( 1.6 ± 0 . 2 ) x 10 3 

(1.5 ± 0.2) x 10 3 

L/S 

0 36"*”° 11 
u.oo-c 09 

0.25 ±0.07 

0-65 ±g ' 24 

0 54 + 0 - 19 
u ' 04t — 0.14 


Note. — The quoted values and uncertainties are 50, 15.87, and 84.13 percentiles of the marginalized MCMC 
posteriors. 

* Value at BJD = 2455633 = 2454833 + 800. 

** Gaussian prior M p /M* = (4.2 ±0.8) x 10 -3 is imposed. The value is based on the average and standard deviation 
of the results by S14, Esteves et al. (2014), and Faigler & Mazeh (2014). 

*** Gaussian prior C = 0.0776 ±0.0200 is imposed. The central value is from the result for n = 3 polytrope by Szabo 
et al. (2012) and the width is chosen to enclose that of the Sun. 


tered on 0.0776 (the value for n = 3 polytrope by Szabo 
et al. 2012) and with the width of 0.02, which is chosen 
to enclose the solar value, 0.059. 

The constraints from the MCMC fit are summarized 
in the middle and bottom parts of Table 2 and the best- 
fit models are plotted with the solid lines in Figure 5. 
Basically, the precession model is compatible with the 
observations both for the light-curve solution and the 
joint solution. The value of J 2 and the corresponding 
precession period, however, are different by a factor of 
a few, in spite of the similar observed slopes in cos?' 0 rb- 
While J 2 = (1.66 ± 0.08) x 10 -4 for the light-curve solu¬ 
tion is consistent with the earlier estimate by Szabo et al. 
(2012), J 2 = (2.1 ± 0.6) x 10 -4 from observed TDVs and 
J 2 = 1.7 x 10 -4 from the stellar model, the joint solution 
yields a smaller value, J 2 = (6.1 ± 0.3) x 10 -5 . 

The difference comes from the different three- 
dimensional architectures of the system described by the 
two solutions. Since all of cosf or bi A, and ?'* are con¬ 
strained from the gravity-darkened light curves, relative 
configuration of the stellar spin and orbital angular mo¬ 
menta are completely specified in three dimensions. This 
means that the phase of the precession during the Kepler 
mission, which corresponds to the left end in the right 
column of Figure 6, is observationally constrained; from 
the top panel, we find that cos ?' 0 rb is closer to the bottom 
of the sine curve for the light-curve solution (blue dashed 
line), while that for the joint solution (red solid line) re¬ 
sides in the phase of a rapid increase. For this reason, 
a larger precession rate (i.e., shorter precession period) 
is required for the light-curve solution to match the ob¬ 
served change in cos/ or b- According to Equations (5) and 
(6), the larger precession rate can be achieved by increas¬ 
ing either J 2 or L/S. However, the larger precession rate 
also induces faster variations in A and i*, contradicting 
their almost constant observed values (middle and bot¬ 


tom panels in Figure 5). The only way to mitigate this 
conflict is to make J 2 larger (i.e., increase the precession 
rate) while keeping L/S small, making it more difficult to 
move stellar spin axis. With Equation (7), this explains 
why Mp/M* is smaller and C is larger for the light-curve 
solution than for the joint solution in Table 2. Accord¬ 
ingly, the bottom panel of the right column in Figure 6 
exhibits the smaller precession amplitude for i* in the 
former solution (blue dashed line) than the latter (red 
solid line). 

The approximately three times difference in the pre¬ 
cession period would be apparent even on the short time 
scale (left column in Figure 6). As shown in the middle 
panel, as large as ~ 10° change in A is expected within 
the next ~ 10 yr for the light-curve solution, which may 
well be detectable given the current precision of the spin 
orbit angle measurement (nominally down to a few de¬ 
grees). On the other hand, A for the joint solution is 
almost constant. From this point of view, the joint solu¬ 
tion may slightly be favored even with the current data, 
because the nearly-constant values observed for A and 
i* are more natural for the joint solution than for the 
light-curve one, for the reasons discussed in the previous 
paragraph. This indication also manifests itself in the 
fact that the resulting M p /M* and C better agree with 
our prior knowledge in the joint solution. 

The more decisive conclusion will be obtained with the 
future follow-up observations of A using Doppler tomog¬ 
raphy, as well as the transit duration observations to bet¬ 
ter constrain cosi or b, and hence the precession rate. If 
our joint solution is correct, variations in A will not be 
detected in near future. On the other hand, if the light- 
curve solution is actually correct and A from the Doppler 
tomography is somehow systematically biased, A should 
change; this temporal variation would be observable with 
the Doppler tomography even if it were biased. Or, it 
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Fig. 5.— Simultaneous fit to the observed cosz or bj A, and z*. (Left) light-curve solution with C 2 = 0. (Right) joint solution. Black points 
are from the light-curve fit (same as Figure 4), and colored solid lines denote the best-fit precession models, which are not the linear fits. 
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year year 

Fig. 6 .— Future evolutions of cosz or b, A, and z* predicted for the best-fit models in Table 2 and Figure 5. From top to bottom, the 
evolutions of cosz or b, A, and z* are plotted for the solution from the light-curve alone (blue dashed line) and joint solution (red solid line) 
obtained with the S14 stellar parameters. The left panels show the short-term (<~ 14 yr, until 2022) behavior, while the right ones are for 
the long-term (~ 1600 yr) variation. 
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may even turn out that the precession scenario is wrong. 
In any case, tracking the future evolution of the system 
configuration can be used for an independent test of our 
solution, not to mention for better constraining stellar 
internal structure via Ji and C, for which few observa¬ 
tional constraints have been obtained. 

5. ANOMALY IN THE TRANSIT LIGHT CURVE OF 

HAT-P-7 

Armed with the methodology established using the dis¬ 
tinct anomaly in Kepler-13A (Section 3), we discuss an¬ 
other, more subtle anomaly in this section. Here the 
methodology is further extended to include the informa¬ 
tion from asteroseismology as well as from the RM effect, 
and applied to an F-type star. 

It has been pointed out in several studies that the 
transit light curve of HAT-P-7 exhibits a small anomaly 
of 0(1O -5 ). Morris et al. (2013), who reported this 
anomaly first, attributed it to the local spot-like grav¬ 
ity darkening induced by the gravity of the Jupiter-mass 
companion HAT-P-7b. They ruled out the gravity dark¬ 
ening of stellar rotational origin on the basis of the in¬ 
spection that the anomaly is localized in a part of the 
transit. Later analyses with more data (e.g., Esteves 
et al. 2013; Van Eylen et al. 2013; Esteves et al. 2014; 
Benomar et al. 2014), however, have shown that the 
anomaly is seemingly correlated over the whole transit 
duration, as in the top panel of Figure 8. Moreover, the 
amplitude of the observed anomaly may be too large to 
be explained by the spot scenario. According to Jackson 
et al. (2012), the planet’s gravity induces the surface tem¬ 
perature variation of “a few 0.1 K,” which leads to the 
surface brightness variation of A F ~ several 100 ppm. 
If a planet crosses over a spot fainter by A F than the 
other part of the stellar disk, amplitude of the expected 
anomaly in the relative flux is about A F x (i? p /i?.*) 2 ~ 
O(ppm), which is order-of-magnitude smaller than the 
observed one. We therefore analyze this anomaly as¬ 
suming that it is originated from the gravity darkening 
induced by stellar rotation, whose effect should not be 
localized but manifest during the whole transit duration. 

Unlike the case of Kepler-13A, anomaly in the tran¬ 
sit light curve is not clear on a quarter-by-quarter ba¬ 
sis for HAT-P-7. In addition, no TTVs/TDVs have 
been detected for this planet. For these reasons, we 
deal with the light curve obtained by folding all the 
available SC, PDCSAP fluxes (QO-17) processed as de¬ 
scribed in Section 2.2. We use the spectroscopic con¬ 
straint usini* = 3.8±1.5kms _1 throughout this section. 
This value is based on Pal et al. (2008), though its er¬ 
ror bar is enlarged to take into account other estimates 
for this quantity that give slightly different values (e.g., 
Winn et al. 2009a). 

5.1. Robustness of the Observed Anomaly 

If the observed anomaly is really due to gravity darken¬ 
ing, it should be persistent over all observation span. It is 
important to confirm the property because Morris et al. 
(2013) only reported the bump before the mid-transit 
time. Thus, we divide the transits into four consecu¬ 
tive groups (Q0-4, 5-9, 10-13, 14 17), phase-fold and fit 
each of them with the model without gravity darkening 
separately, and examine the shapes of the residuals. Al¬ 


though fewer numbers of transits lead to noisier phase- 
folded light curves, ten-minute binned residuals in the 
left column of Figure 7 exhibit a similar feature (bright¬ 
ening before mid-transit and dimming after it) in every 
span of data. 

Besides, Van Eylen et al. (2013) reported seasonal vari¬ 
ation in the transit depth depending on the quarter, 
which is reproduced in our analysis with QO-17 data. 10 
To confirm that the anomaly is not an artifact related to 
this seasonal variation, we also perform a similar analysis 
as above but this time grouping the transits that have 
similar depths. As shown in the right column of Figure 
7, we find that the same feature is apparent regardless 
of the season and the anomaly is not affected by the sys¬ 
tematic depth variation. For this reason, along with its 
unconstrained origin, we do not try to make corrections 
for this systematic in the following analyses. 

5.2. Results 

As in Section 3, we consider both light-curve solution 
and joint solution that takes into account the constraints 
from other observations. First, the light-curve solution 
is obtained with c-i fixed to be zero (Figure 8, second 
and third columns in Table 3). In this case, we find 
two solutions with different signs of cosi 0 rb> which are 
indistinguishable in terms of the minimum x 2 . 11 The 
values quoted in Table 3 are the median, 15.87, and 84.13 
percentiles of the MCMC posteriors sampled with emcee 
(Foreman-Mackey et al. 2013). 12 Our model reasonably 
reproduces the global feature of the anomaly (positive 
before the mid transit and negative after it), yielding 
Ay 2 — 166 for ~ 420 degrees of freedom. We compute 
the Bayesian information criterion (BIC) for the best- 
fit models with and without gravity darkening, and find 
ABIC = 129, which formally indicates that the gravity- 
darkened model is strongly favored. 

Our solution points to a nearly pole-on configuration 
with i* ~ 0°. This conclusion is consistent with the 
recent asteroseismic analyses by Benomar et al. (2014) 
and Lund et al. (2014), but the nominal constraint on t*. 
from the gravity-darkened model is much tighter. On the 
other hand, A is not constrained very well with the light 
curve asymmetry alone. The difficulty is inevitable in the 
pole-on configuration, where the brightness distribution 
on the stellar disk is almost axisynnnetric even in the 
presence of gravity darkening. In such a case, ip is always 
close to 90° regardless of A. 

One remaining issue regarding our solution is that the 
resulting rotation frequency may be too large. Given 
the age (~ 2Gyr) and B — V (= 0.495 ± 0.022; Lund 

10 We also reported a similar phenomenon in Kepler-13A; see 
Section 4.1 and Figure 4. 

11 The existence of the two solutions in this case should be dis¬ 
tinguished from the degeneracy intrinsic to the gravity-darkening 
method. For each of the two solution listed here, there addition¬ 
ally exists the model that yields exactly the same light curve, 
where cos i or n is replaced with — cos i or b and A with it — A. These 
intrinsically-degenerate solutions are not discussed here because 
they are in any case rejected in the joint solution, where A is con¬ 
strained by the prior. This is the same logic as used in the last 
paragraph of Section 3.1. 

12 We also applied the residual permutation method described 
in Winn et al. (2009b) for another estimate of the parameter un¬ 
certainties, and confirmed that they are not significantly affected 
by the correlated noise component. 
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et al. 2014) of the host star, the rotation frequency from 
the light-curve solution, / rot = 7.7 ± 0.2 yMlz (equivalent 
to P ro t ~ 1.5days), is consistent with the gyroclrronol- 
ogy relation by Meibom et al. (2009); see Section 6 of 
Lund et al. (2014). However, our value of f m t is much 
larger than those from asteroseismology, OTO/g/g ^zHz 
(68% credible interval by Benomar et al. 2014) and 
< 0.8748/iHz (lcr upper limit by Lund et al. 2014). In 
fact, the prior used in these analyses, |/ ro t| < 8/xHz, 
does not fully cover the range we investigate here with 
the gravity-darkened light curve. Still, the discrepancy is 
only weakly reduced even with the new analysis adopt¬ 
ing the prior range extended up to 17 yHJz, which yields 
/rot = 0.82/q;°q ^zHz as the 68% credible interval (by 
courtesy of Othrnan Benomar; see also Benomar et al. 
2014). 

To examine if the gravity-darkened model is compati¬ 
ble with the seismic analysis, we then search for a joint 
solution including the constraints both from the RM 
measurement and asteroseismology. From the RM ef¬ 
fect, we incorporate the constraint A = 172° ± 32°, 
which comes from the average and standard deviation 
of the analyses for the three different radial velocity data 
(Benomar et al. 2014). From asteroseismology, we adopt 
the above updated posterior for / rot as the prior, and 
performed an MCMC sampling with emcee. To prop¬ 
erly take into account the uncertainty from the limb- 
darkening profile, C 2 is also floated. The resulting cred¬ 
ible intervals are summarized in the fourth and fifth 
columns in Table 3, and the model that maximizes the 
likelihood multiplied by the prior on / rot is plotted with 
a dashed line in Figure 8. We again find two equally 
good solutions, both of which indicate nearly pole-on 
configurations with slightly prograde and retrograde or¬ 
bits, ip = 101° ± 2° and ip = 87° ± 2°. Although the 
resulting / rot still prefers a higher rotation rate than that 
from asteroseismology, their difference is now mitigated 
to the 2 <t level; we construct the probability distribution 
for A/ rot , / r ot from out joint analysis minus f vo t from 
asteroseismology, using their posteriors and find its 2er 
credible region as A/ rot = 4 . 9 / 5 ;° l^Hz. We argue that 
the level of discrepancy is acceptable, considering that 
the rotational mode splitting is not clearly detected in 
the power spectrum of HAT-P-7’s light curves. 

Finally, it is also worth considering the case where 
[3 ^ 0.25, given the unconstrained nature of the grav¬ 
ity darkening in F dwarfs. Smaller values of (3 ~ 0.08 
are usually expected for solar-like stars with convective 
envelopes (e.g., Lucy 1967; Claret 1998), while Espinosa 
Lara & Rieutord (2011) and Espinosa Lara & Rieutord 
( 2012 ) argue that /3 is close to 0.25 in the limit of slow 
rotation under several assumptions. We repeat the above 
joint analysis floating (3 with the prior uniform between 
0 and 0.3, and obtain (3 = 0. 26/^05 for both solutions 
in Table 3. On the one hand, the fact may support the 
claims by Espinosa Lara & Rieutord (2011) and Espinosa 
Lara & Rieutord (2012); on the other hand, it may simply 
indicate some incompleteness in our gravity-darkening 
model, as also suggested by the tension in / rot and the 
still correlated residuals before the mid transit (bottom 
panel of Figure 8 ). Indeed, if (3 = 0.08 is adopted, we 
find that even higher rotation rate (> 10 /zHz) is favored, 
making the discrepancy with asteroseismology more se¬ 


rious. Although the validity of /3 we obtain is beyond 
the scope of this paper, we note that our conclusion for a 
pole-on orbit is robust against the adopted value of /3; in 
both analyses where /? is fitted and f3 is fixed to be 0.08, 
the constraints on ip differ less than lcr from the results 
in Table 3. 

6. SUMMARY 

6.1. Kepler-13A 

First, we analyze the gravity-darkened transit light 
curve of Kepler-13A adopting the same model and stellar 
parameters as in the previous study by Bll. We repro¬ 
duce the spin-orbit angles obtained by Bll with more 
data (called “light-curve solution” in this paper) and also 
find that the choice of the stellar mass, stellar effective 
temperature, usinz*, or contaminated flux affects A or 
z* by less than about 10°. If we fit C 2 = u\ — 112 as well 
as ci = Ui + U 2 in the quadratic limb-darkening law, on 
the other hand, a broader range of the spin-orbit an¬ 
gle is allowed. In fact, this additional degree of freedom 
may explain the discrepancy between the solution by Bll 
and the Doppler tomography result by Johnson et al. 
(2014). Our new “joint solution” includes z* = 81° ± 5°, 
A = -59° ± 2°, ip = 60° ± 2°, and P rot = 24 ± 2hr. 
Although the joint solution is compatible with all of the 
observations made so far, introducing additional free pa¬ 
rameter C 2 is not statistically justified, nor is it clear if 
the best-fit value for C 2 is physically plausible. 

To examine the above issues from a dynamical point 
of view, we also analyze the spin-orbit precession in this 
system. By analyzing the light curves from each quar¬ 
ter separately, we confirm that the variation in | cos z or b | 
causes the transit duration variations first reported by 
Szabo et al. (2012), with more elaborate model taking 
into account the gravity darkening. This variation is con¬ 
sistent with the precession of the stellar spin and orbital 
angular momenta around the total angular momentum 
of the system, induced by the oblateness of the rapidly 
rotating host star. We thus fit each transit with the 
gravity-darkened model to determine cosz 0 rb, A, and z* 
as a function of time, and then fit them with the pre¬ 
cession model to constrain the stellar quadrupole mo¬ 
ment J 2 . For the light-curve solution and the joint solu¬ 
tion, we respectively find J 2 = (1.66 ± 0.08) x 10 -4 and 
J 2 = (6.1±0.3) x 10~ 5 , which are different by a factor of a 
few. Our results predict detectable variations in A on 10- 
yr timescale for the light-curve solution, while it should 
be almost constant for the joint solution. The difference 
suggests that the future follow-up observations can be 
used to confirm or refute the joint solution we proposed, 
as well as to improve the constraint on J 2 . 

6.2. HAT-P-7 

Although the anomaly in the transit light curve is much 
more subtle compared to Kepler-13Ab, we confirm that 
the asymmetric residual (not only the bump reported 
by Morris et al. (2013) but also the dip) exists continu¬ 
ously in the transits of HAT-P-7b. Thus, we perform the 
analysis assuming that the gravity-darkening is a viable 
explanation for the anomaly. Gravity-darkened transit 
model favors a nearly pole-on orbit (ip = 101° ± 2° or 
ip = 87° ±2°) and the gravity-darkening exponent [3 close 
to 0.25. The constraint on ip is insensitive to the choice of 
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Fig. 7.— Robustness of the detected anomaly. Residuals of the model fits (without gravity darkening) to the phase-folded transit light 
curves are plotted. Gray dots are residuals for the one-minute binned data, and black ones with error bars are the residuals averaged into 
ten-minutes bins. Vertical dashed and dotted lines correspond to the beginnings and ends of the ingress and egress. (Left column) Transits 
folded over different epochs. From top to bottom, light curves from Quarters 0—4, 5—9, 10—13, and 14-17 are folded. (Right column) 
Transits grouped by the CCD module used to observe the target. From top to bottom, light curves taken with the modules 17, 19, 9, and 
7 are folded. 


the limb-darkening parameters or the gravity-darkening 
exponent. 

On the other hand, the stellar rotation rate from the 
gravity-darkening analysis is about 2cr higher than the 
value from asteroseismology. In addition, the value of 
/3 ~ 0.25 we obtained may be too large for a star with a 
convective envelope. These facts, as well as the subtle¬ 
ness of the detected anomaly, may suggest some incom¬ 
pleteness in the current modeling or other origins for the 
anomaly, and should be addressed in future studies. 

7. CONCLUSION 

Our present analysis reproduces the results by Bll 
with more data and thus strengthens the reliability of the 
gravity-darkening method for the spin-orbit angle deter¬ 
mination. In contrast, we also find that the spin-orbit 
angle obtained from the gravity-darkened transit light 
curve strongly depends on the assumed limb-darkening 
profile. Depending on its choice, the resulting spin-orbit 
angle can vary by several tens of degrees. Thus, the reli¬ 
able modeling of the limb-darkening effect is crucial for 
this method. 

Nevertheless, if A is constrained from other observa¬ 
tions, i* is well determined along with the limb-darkening 
parameters. Hence the gravity-darkening method still 
provides valuable information on the true spin-orbit an¬ 


gle tp, which is complementary to A from the RM effect or 
Doppler tomography. Indeed, such an example is already 
seen in an eclipsing binary system DI Her (Philippov & 
Rafikov 2013). In addition, synergy with asteroseismol¬ 
ogy is also promising because it constraints / ro t and '<*, 
which are both essential in the modeling of gravity dark¬ 
ening. The joint analyses of these kinds may in turn help 
us to better understand the mechanisms of gravity dark¬ 
ening itself, since they enable the measurements of /3 for 
stars not in close binary systems and hence free from the 
strong tidal distortion. 

If combined with continuous, high-precision photom¬ 
etry as achievable with space-borne instruments, the 
gravity-darkening method also provides a way to mon¬ 
itor the angular momentum evolution in the system. 
Modeling of the spin-orbit precession allows us to ac¬ 
cess the internal structure of the rotating star through 
its quadrupole moment or moment of inertia. It is also 
possible to determine the three-dimensional configura¬ 
tion of the system from a dynamical point of view (c.f., 
Philippov & Rafikov 2013; Barnes et al. 2013). Such in¬ 
formation will be valuable in simulating the dynamical 
histories of individual systems to decipher the origin of 
the spin-orbit misalignment. 

The author is grateful to the entire Kepler team for 
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Fig. 8. — Fitting the gravity-darkened model to the phase-folded transit of HAT-P-7b. The meanings of the symbols are the same as those 
in Figure 2, but this time the joint solution incorporates the constraints on A from the RM measurement and on / ro t from asteroseismology. 
The light-curve solution and joint solution are almost indistinguishable in this case, as expected from the similar values of \ 2 . 
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TABLE 3 

Results for the transit of HAT-P-7b 


light-curve solution (C 2 = 0) joint solution (02 fitted) 

Solution 1 Solution 2 Solution 1 Solution 2 


( Constraints) 
usini* (kins -1 ) 

A 

(Fitted Parameters) 

Mi, (M 0 ) 

T*, p ole (K) 

p* (gem 3 ) 

ci 

C2 

t c (10 -5 day)* 

P (day) 

COS ^orb 

Rp/R* 

Fo 

/rot (mHz) 

, V*** 

U (deg) 

A (deg) 

0 

(Derived Parameters) 
P ro t (day) 
i> ( de g) 

impact parameter 
oblateness 


3.8 ± 1.5 


3.8 ± 1.5 


3.8 ± 1.5 
172 ± 32 


3.8 ± 1.5 
172 ± 32 


1.59 (fixed) 
6310 (fixed) 
0.2789 ± 0.0006 
0.498 ± 0.003 
0 (fixed) 
-1.5 ±0.4 


-0.1195 ±0.0004 
u.u< / o/_ 0 00009 


7.7 ±0.2 
3.3j 

133 - w 
0.25 (fixed) 


5 + 1.2 
3 - 1.0 
o+19 


1.59 (fixed) 
6310 (fixed) 
0.2789 ± 0.0006 
0.498 ± 0.003 
0 (fixed) 
-1.5 ±0.4 


1.59 (fixed) 
6310 (fixed) 
0.2790 ± 0.0005 

0-507l°;°?6 

0.07l°;°2 

-1.6 ±0.4 


■ 2.204735471 (fixed) 


0.1195 ±0.0004 -0.1194 ±0.0003 

0.07759 ± 0.00003 
0.9999998 ± 0.0000005 - 


0 077^7+®'00005 
u.u/ _ Q 00009 


7.7 ±0.2 
o q+1.3 

49^21 
0.25 (fixed) 


+ 2 . 6 ** 
1.7 
+3.3 
- 2.0 
+ 12 


6.1 

5.3 

142_i6 
0.25 (fixed) 


1.59 (fixed) 

6310 (fixed) 

0.2784 ± 0.0005 

0 508"*” 0 ' 007 

u.ouo_ 0 015 

0 - 08 io.il 

-1.1 ±0.4 


0.1198 ±0.0003 
0 07749~*~°' 00003 

u.ui 00004 

r £. + 2.4** 
5.6_1 7 
r 0 + 3.7 
0 . 0-21 

136122 

0.25 (fixed) 


1.51 ±0.03 
99^4 

0.496 ± 0.001 


1.51 ±0.03 
81^2 

0.496 ± 0.001 


0.0149 ± 0.0006 0.0149 ± 0.0007 


1.9 


+0.7 

- 0.6 


101 ± 2 

0.496 ± 0.001 

0 009 +0 01 ° 
u.uuj_ 0 005 


2.1 


+0.9 

- 0.6 


87 ± 2 

0.497 ± 0.001 

0.008ig;“t 


Xmin/ dof 453/424 455/424 450/424 451/424 


Note. — The quoted values and uncertainties are 50, 15.87, and 84.13 percentiles of 
the marginalized MCMC posteriors. For the light-curve solution, Xmin is the value of + 
computed from Equation (4) for the maximum likelihood model. Equation (4) is also used 
for the joint solution, but Xmin in this case is computed for the model that maximizes the 
likelihood function multiplied by the prior on / ro t- 

* Measured from the transit epoch fo(BJD — 2454833) = 120.358522 ± 0.000005 obtained 
with the transit model without gravity darkening. 

** Posterior from the seismic analysis is used as the prior. 

*** We impose the prior uniform in cosi*, rather than in i*, which corresponds to the isotropic 
distribution for the spin direction. 
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